function [simdout] = simdata(vardata,...
    nlags,nstep,ndraws,trimpzout,trerzout,trazero,tra0betaz,exo);

uz = trerzout;
betavar = tra0betaz';
a0inv = inv(trazero);
[nobs,nvars] = size(vardata); 
sizeimp = [size(trimpzout)  ndraws];
nobs = nobs-nlags;
nexo=size(exo,2);

%do simulations. 
 for icnt =1:ndraws;
   %if rem(icnt,ndraws/100*50)==0;
   %    disp([int2str(icnt/ndraws*100) '% of bootstrapped data generated...']);
   %end
   erz = uz(ceil(size(uz,1)*rand(nobs,1)),:);
   %I generate a random number bounded between 0 and # of residuals
   %I then use the ceil function to select that row of the residuals
   %this is equivalent to sampling with replacement.    
   ydata = zeros(nobs+nlags,nvars);
   ydata(1:nlags,1:nvars) = vardata(1:nlags,1:end); %use the original values as starting values
   
   %set up lagged values. 
   zlag = reshape(ydata(nlags:-1:1,1:nvars)',1,nlags*nvars);
   %generate simulated data. 
   for jj=nlags+(1:nobs);
      if size(exo,2)>0;
        ydata(jj,:)=exo(jj,:)*betavar(1:nexo,:)+zlag*betavar(nexo+1:size(betavar,1),:)+(a0inv*erz(jj-nlags,:)')';
      else
        ydata(jj,:)=zlag*betavar(1:size(betavar,1),:)+(a0inv*erz(jj-nlags,:)')';
      end
        zlag = [ydata(jj,:) zlag(1,1:(nlags-1)*nvars)];
   end;
   simdout(:,:,icnt) = ydata;
 end